/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		material_interface.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-02-18 13:56
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-07-20 13:00
*
*   Description:
*
=============================================================================*/
#if !defined(INTERFACE_FOR_MATERIAL_CLASS_H)
#define INTERFACE_FOR_MATERIAL_CLASS_H

#include "../arithmetic/tensor.hpp"
#include "../types.hpp"
#include "../util/factory.hpp"
#include "../util/indexed.hpp"
#include "../util/record.hpp"
#include <boost/property_tree/ptree.hpp>

namespace geoxfem {
namespace material {

template <typename DimensionType>
struct get_dee_matrix { /*{{{*/
  static MatrixXR apply(const _real&, const _real&) {
    GEOXFEM_NOT_IMPLEMENTED();
  }
};

template <>
MatrixXR get_dee_matrix<geometry::plane_strain>::apply(const _real& E,
                                                       const _real& v) {
  MatrixXR dee;
  dee.setZero(3, 3);

  _real k0  = v / (r1 - v);
  _real G   = (r1 - r2 * v) / (r1 - v) * r_5;
  _real fac = E * (r1 - v) / (r1 + v) / (r1 - r2 * v);

  dee << r1, k0, r0,  //
    k0, r1, r0,       //
    r0, r0, G;

  return dee * fac;
}
template <>
MatrixXR get_dee_matrix<geometry::plane_stress>::apply(const _real& E,
                                                       const _real& v) {
  MatrixXR dee;
  dee.setZero(3, 3);

  _real k0  = v;
  _real G   = (r1 - v) * r_5;
  _real fac = E / (r1 - v * v);

  dee << r1, k0, r0,  //
    k0, r1, r0,       //
    r0, r0, G;

  return dee * fac;
}
template <>
MatrixXR get_dee_matrix<geometry::three_dimension>::apply(const _real& E,
                                                          const _real& v) {
  MatrixXR dee;
  dee.setZero(6, 6);

  _real k0  = v;
  _real G   = (r1 - v) * r_5;
  _real fac = E / (r1 - v * v);
  dee << r1, k0, r0, r0, r0, r0,  //
    k0, r1, r0, r0, r0, r0,       //
    r0, r0, r1, r0, r0, r0,       //
    r0, r0, r0, G, r0, r0,        //
    r0, r0, r0, r0, G, r0,        //
    r0, r0, r0, r0, r0, G;

  return dee * fac;
}

/*}}}*/

template <typename DimensionType>
struct MaterialInterface : util::Indexed {
  using Material = MaterialInterface<DimensionType>;
  using ptree    = boost::property_tree::ptree;
  using Record   = util::iRecord<DimensionType>;
  using Stress   = geoxfem::arithmetic::Stress<DimensionType>;
  using Strain   = geoxfem::arithmetic::Strain<DimensionType>;
  using factory  = dp14::factory<Material, const ptree&>;
  static constexpr size_t dimension = traits::dimension<DimensionType>::value;

  /*constructor{{{*/

  /// default constructor
  MaterialInterface(const ptree& data) {
    m_param.addField("volume_weight", data.get("volume_weight", _real{19}));
    m_param.addField("friction_angle", data.get("friction_angle", _real{30}));
    m_param.addField("dilation_angle", data.get("dilation_angle", _real{30}));
    m_param.addField("cohesion", data.get("cohesion", _real{10}));
    m_param.addField("zeta", data.get("zeta", _real{0.1}));
    m_param.addField("K_IC", data.get("K_IC", _real{2.9}));
    m_param.addField("elastic_modulus",
                     data.get("elastic_modulus", _real{10000}));
    m_param.addField("elastic_poisson",
                     data.get("elastic_poisson", _real{0.3}));

    Stress perm;
    for(Eigen::Index i = 0; i != 3; ++i) {
      std::string name = "perm_x" + std::to_string(i);
      perm(i, i)       = data.get(name, _real{1e-5});
    }
    m_param.template addField<Stress>("perm", perm);
  }

  /**
   * @brief ~MaterialInterface
   */
  virtual ~MaterialInterface() = default;
  /*}}}*/
  /*misc{{{*/

  static Material::factory& getFactory() {
    return Material::factory::instance();
  }
  /// get class name
  virtual std::string getClassName() const {
    return "MaterialInterface";
  }

  /**
   * @brief operator<<
   *
   * @param os
   * @param mat
   *
   * @return
   */
  friend std::ostream& operator<<(std::ostream& os, const Material& mat) {
    os << "\nMaterial (" << mat.m_type << ")\n"
       << "-------------\n";
    os << mat.m_param;
    os << "-------------\n";
    return os;
  }
  /*}}}*/
  /*member manuvar{{{*/

  /**
   * @brief getParams
   *
   * @return
   */
  Record getParams() const {
    return m_param;
  }
  /**
   * @brief &getParams
   *
   * @return
   */
  Record& getParams() {
    return m_param;
  }
  /**
   * @brief getParam
   *
   * @param name
   *
   * @return
   */
  template <typename vtype = std::string>
  vtype getParam(const std::string& name) const {
    return this->m_param.template get<vtype>(name);
  }
  /*}}}*/
  /* matrix manuvar {{{*/
  /** get seepage matrix according to perm_x, perm_y and perm_z
   *
   * @return
   */
  MatrixXR getSeepageTensor() const {
    const auto& perm = m_param.template getValue<Stress>("perm");
    return perm.topLeftCorner(dimension, dimension);
  }

  /** get [D_{elastic}] matrix according to given E and v
   *
   * @param E
   * @param v
   * @return
   */
  MatrixXR getElasticDeeMatrix(const _real& E, const _real& v) const {
    return get_dee_matrix<DimensionType>::apply(E, v);
  }
  //________________________________________________________________

  /** get [D_{elastic}] matrix according to Linear E and v
   *
   * @return
   */
  MatrixXR getElasticDeeMatrix() const {
    _real E = m_param.template getValue<_real>("elastic_modulus");
    _real v = m_param.template getValue<_real>("elastic_poisson");
    return getElasticDeeMatrix(E, v);
  }
  /*}}}*/
  /*non-linear realted{{{*/

  /** integrate sigma through epsilon path with substepping method
   *
   * @param de  :
   * @param s                :
   * @param param
   */
  void
  IntegralThroughEpsilonPath(const Strain& de, Stress& s, Record& param) const {
    /* Reference: */
    /* @article{ */
    /* doi:10.1002/nme.1620240505, */
    /* author = {Sloan, S. W.}, */
    /* title = {Substepping schemes for the numerical integration of
     * elastoplastic stress–strain relations}, journal = {International Journal
     * for Numerical Methods in Engineering}, volume = {24}, number = {5}, pages
     * = {893-911}, doi = {10.1002/nme.1620240505}, url =
     * {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620240505}, eprint
     * = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620240505}, */
    /* abstract = {Abstract This paper describes two substepping schemes for
     * integrating elastoplastic stress–strain relations. The schemes are
     * designed for use in finite element plasticity calculations and solve for
     * the stress increments assuming that the strain increments are known.
     * Both methods are applicable to a general type of constitutive law and
     * control the error in the integration process by adjusting the size of
     * each substep automatically. The first method is based on the well-known
     * modified Euler scheme, whereas the second technique employs a high order
     * Runge–Kutta formula. The procedures outlined do not require any form of
     * stress correction to prevent drift from the yield surface. Their utility
     * is illustrated by analysis of typical boundary value problems.}, */
    /* year = {1987}  */
    /* } */

    _real precision = _real(1e-2);

    /* 1. trial delta sigma */
    Stress ds;
    ds.fromVoigt(getElasticDeeMatrix() * de.Voigt());

    /* 2. if trial sigma lies in elastic zone */
    if(getYieldFunction(s + ds, param) <= 0)
      s += ds;
    /* 3. if trial sigma lies in plastic zone */
    else {
      /* 2.1 get alpha: the ratio of elastic/all in trial delta sigma */
      const _real alpha = getRatioOfNonLinear(s, ds, param);
      /* dpe: delta plastic epsilon */
      const Strain dpe = (1 - alpha) * de;
      /* 2.2 add delta elastic sigma.
       * find the elastic sigma at yield surface.  */
      s += ds * alpha;
      /* 2.3 iteration to find final sigma/epsilon */
      for(_real T = 0, dT = 1; T < 1;) {
        dT = (dT < 1 - T) ? dT : 1 - T;
        /* ddpe: trial detal plastic epsilon */
        Strain ddpe = dT * dpe;
        /* 2.3.1. first trial */
        Stress ds1;
        ds1.fromVoigt(getDeeMatrix(s, param) * ddpe.Voigt());
        /* 2.3.2. get new sigma and state paramenter
         * (refers to new yield function) */
        Stress s1     = s + ds1;
        Record param1 = param;
        updateStateParemeter(s1, ddpe, param1);
        /// 2.3.3. second trial
        Stress ds2;
        ds2.fromVoigt(getDeeMatrix(s1, param1) * ddpe.Voigt());
        /// 2.3.4. average
        Stress ds_mid = (ds1 + ds2) * 0.5;
        /// 2.3.5. convergence check
        Stress diff          = (ds2 - ds1) * 0.5;
        _real  relative_diff = diff.norm() / (s + ds_mid).norm();
        if(relative_diff < precision) {
          T += dT;
          /// add plastic delta sigma
          s += ds_mid;
          updateStateParemeter(s, ddpe, param);
        }
        else {
          _real beta = std::sqrt(precision / relative_diff) * _real(0.8);
          beta       = beta > _real(0.1) ? beta : _real(0.1);
          beta       = beta < _real(2.0) ? beta : _real(2.0);
          dT *= beta;
        }
      }
    }

    return;
  }

  _real getRatioOfNonLinear(const Stress& s,
                            const Stress& ds,
                            const Record& param,
                            const _real&  precision = _real(0.01)) const {
    /// Reference:
    /// @ misc{ports2001finite,
    /// title={Finite element analysis in geotechnical engineering},
    /// author={Ports, DM and Zdravkovic, L},
    /// year={2001},
    /// publisher={Thomas Telford, London}
    /// }
    _real _alpha_last = 0;
    _real _alpha_curr = 1;
    _real _alpha_next;

    /// secant iteration
    while(std::abs(_alpha_last - _alpha_curr) > precision) {
      _real _f_curr = getYieldFunction(s + _alpha_curr * ds, param);
      _real _f_last = getYieldFunction(s + _alpha_last * ds, param);

      if(_f_last == _f_curr)
        _alpha_next = 0;
      else
        _alpha_next =
          _alpha_curr
          - _f_curr / (_f_curr - _f_last) * (_alpha_curr - _alpha_last);

      _alpha_last = _alpha_curr;
      _alpha_curr = _alpha_next;
    }

    return _alpha_curr;
  }

  /**   get non-linear dee matrxi according to given element state parameter
   * refers to [elastic dee matrix] in elastic material or [non-linear elastic
   * dee matrix] in hypo-elastic material or [elasto-plastic dee matrix] in
   * elasto-plastic material1
   *
   * @param s                :
   * @param param
   * @return
   */
  virtual MatrixXR getDeeMatrix(const Stress& s, const Record& param) const = 0;

protected:
  /** update stress history related parameters for given sigma and epsilon
   * (only override in hardening/softening material for changing yield
   * function)
   *
   * @param s                 :
   * @param dpe
   * @param param        :
   */
  virtual void updateStateParemeter(const Stress& s,
                                    const Strain& dpe,
                                    Record&       param) const = 0;
  //________________________________________________________________

  /** get the value of the yield function
   *
   * @param s                :
   * @param param
   * @return
   */
  virtual _real getYieldFunction(const Stress& s,
                                 const Record& param) const = 0;
  /*}}}*/

protected:
  Record      m_param;
  std::string m_type;
};

}  // namespace material
namespace traits {

template <typename GaussOrCell>
struct get_material {
  static const material::MaterialInterface<typename GaussOrCell::DimensionType>*
  apply(const GaussOrCell& obj) {
    return obj.getMaterial();
  }
};

}  // namespace traits
}  // namespace geoxfem

#endif  // _INTERFACE_FOR_MATERIAL_CLASS_H
